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§^ ; Abstract 

p • We compare the leading and next-to- leading order QCD predictions for the flux of 

atmospheric muons and neutrinos from decays of charmed particles. We find that the 
D ". full NLO lepton fluxes can be approximated to within ~ 10% by the Born-level fluxes 
\ multiplied by an overall factor of 2.2 — 2.4, which depends slightly on the PDF. This 
supports the approach in Thunman, Ingelman, Gondolo (1996). We also find that their 
^ very low lepton fluxes are due to the mild slope they used for the gluon distribution 
. function at small momentum fractions, and that substantially larger lepton fluxes result 
when the slope of the gluon distribution function at small momentum fractions is 
larger. 



1 Introduction 



The flux of atmosplieric neutrinos and muons at very liigli energies, above 1 TeV, 
passes from being originated in the decays of pions and kaons to being predominantly 
generated in semileptonic decays of charmed particles (see for example [Q). This flux is 
of importance for large area detectors of high energy cosmic neutrinos. Future km^ ar- 
rays would be able to observe muons and neutrinos with energies that may reach 10^^ 
GeV. Atmospheric muons and neutrinos would be one of the most important back- 
grounds, limiting the sensitivity of any "neutrino telescope" to astrophysical signals. 
Besides, they might be used for detector calibration and perhaps, more interestingly, 
be exploited to do physics, e.g. study neutrino masses. 

Present experimental attempts to detect atmospheric muons from charm are spoiled 
by systematic errors. Theoretical predictions depend strongly on the reliability of the 
model adopted for charm production and decay and differ by orders of magnitude, due 
to the necessity of extrapolating present accelerator data on open charm production 
in fixed target experiments, at laboratory energies of about 200 GeV, to the larger 
energies needed for atmospheric neutrinos, from 10'^ to 10^ GeV (at about 10^ GeV the 
rates become too small for a km^ detector). These energies, from 40 GeV to 14 TeV in 
the center of mass, are comparable to the energies of the future RHIC at Brookhaven, 
200 GeV, and LHC at CERN, 7 TeV. 

The theoretically preferred model, perturbative QCD (pQCD), was thought to be 
inadequate because it could not account for several aspects of some of the early data 
on open charm production (in conflict with each other, on the other hand [Q), and 
because of a sensitivity of the leading-order (LO) calculation, the only existing until 
recently, to the charm quark mass, to the low partonic momentum fraction, x, behavior 
of the parton distributions and to higher order corrections. So, even if some now- 
obsolete pQCD calculations have appeared H, the models for charm production 
traditionally favored in studies of atmospheric fluxes have been non-perturbative: for 
example, besides semi-empirical parametrizations of the cross section, the quark-gluon 
string model (QGSM, a.k.a. dual parton model), based on Regge asymptotics, and 
the recombination quark-parton model (RQPM), incorporating the assumption of an 
intrinsic charm component in the nucleon (see 0). 

Today, however, pQCD predictions and experimental data are known to be compat- 
ible 0, H, ^,0: charm production experiments form a consistent set of data, and the 
inclusion of next-to-leading order (NLO) terms has been a major improvement over the 
leading-order treatment. Quoting from Appel 0, "the success of these calculations has 
removed the impetus to look for unconventional sources of charm production beyond 
the basic QCD" . 

A study based on pQCD was therefore performed in Ref. |^ (called TIG from 
now on). CLEO and HERA results were incorporated, but for simplicity the LO 
charm production cross section was adopted, multiplied by a constant K factor of 2 
to bring it in line with the next-to-leading order values, and supplemented by parton 
shower evolution and hadronization according to the Lund model. The neutrino and 
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muon fluxes from cliarm were found to be lower than the lowest previous prediction, 
namely a factor of 20 below the RQPM pl, of 5 below the QGSM |T|, 0], and of 3 
below the lowest curve in Ref. IQ. 

Here we use the same treatment of TIG, except for the very important difference 
of using the actual next-to-leading order pQCD calculations of Mangano, Nason and 
Ridolfi |T^ (called MNR from now on), as contained in the program we obtained 
from them (see also [0), to compute the charm production cross sections. These are 
the same calculations used currently to compare pQCD predictions with experimental 
data in accelerator experiments. The main goal of this paper is to compare the fluxes 
obtained with the NLO and with the LO, i.e. we will compute the K factor for the 
neutrino and muon fluxes. This K factor is necessarily different from the K factor 
for charm production (which can be found in the literature), because only the forward 
going leptons contribute significantly to the atmospheric fluxes. 

A similar comparison was very recently made in [|17| , using the approximate analyt- 
ical solutions introduced by TIG to the cascade equations in the atmosphere. We make 
instead a full simulation of the cascades, using the combined MNR and PYTHIA pro- 
grams. These two treatments of the problem are complementary. For comparison, we 



include results obtained with the CTEQ 3M gluon structure function used in Ref. [17 



We find our CTEQ 3M results to be close to those of the PRS study, in spite of the 
very different approaches used in the two calculations. 

Addressing right away a concern that has been expressed to us several times, about 
the applicability of perturbative QCD calculations, mostly done for accelerator physics, 
to the different kinematic domain of cosmic rays, we would like to point out that, since 
the characteristic charm momentum in our simulations is of the order of the charm 
mass, k ~ 0(mc), we do not have here the uncertainty present in the differential cross 
sections [|T3|, when is much larger than rric (as is the case in accelerators), due to 
the presence of large logarithms of (fcf^ -|- m^)/m^. Depending on the steepness of the 
gluon structure function we take, we do have, however, large logarithms, known as 
"ln(l/x)" terms, where x ~ ^jAml/ s (s is the hadronic center of mass energy squared) 
is the average value of the hadron energy fraction needed to produce the cc pair. These 
should not be important for steep enough gluon structure functions (namely for values 
of A in Eq. (9) not very close to zero), but we have not made any attempt to deal with 
this issue. 

In the next section of this paper we explain our normalization of the NLO charm 
production cross section in the MNR program. In Sect. |^ we describe the computer 
simulations used to calculate the neutrino and muon fluxes. In Sect. ^ we show the 
results of our simulations, we discuss the differences between a NLO and a LO approach 
and we make a comparison with the fluxes of the TIG model. 

In this paper we consider only vertical showers for simplicity (the same was done 
by TIG). 
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2 Charm production in perturbative QCD 



In this section, we show evidence that perturbative QCD gives a fair description of 
the present accelerator data on open charm production in the kinematic region most 
important for cosmic ray colhsions in the atmosphere. 

There are still not many experiments on open charm production with good enough 
statistics, despite the recent improvements, but many are expected in the near future. 

We use a NLO approach which is based on the MNR calculation, for which we have 
obtained the computer code. The NLO cross section for charm production depends on 
the choice of the parton distribution functions (PDFs) and on three parameters: the 
charm quark mass rric, the renormalization scale n^, and the factorization scale hf- 



2.1 Choice of rric, Hr, Hf 

MNR have two default choices of rric, fJ^R and fip- for total cross sections they choose 
rric = 1.5 GeV, = rric, fJ^F = '^ruc, for differential cross sections they choose instead 



TTic = 1.5 GeV, fiR = rriT, ^if = '^rriT, where m-r = yk"^ + rnl is the transverse mass 
The current procedure to reproduce the measured differential cross sections [^, ^, |10 



is to use the MNR default choices for these three parameters and multiply the result 
by the global factor of about 2 or 3 necessary to match the predicted and measured 
total inclusive cross sections. Although this procedure might be acceptable in face of 
the uncertainties in the pQCD predictions, we find it unsatisfactory from a theoretical 
point of view. We prefer to fit the differential and total cross sections with one and 
the same combination of rric, fJ'R, and fip- 

We make separate fits of rric, fJ'R, and fip for each of the following sets of PDFs: 



MRS Rl, MRS R2 p|, CTEQ 3M [|l|] and CTEQ 4M (see the next subsection 
for details). 

We are aware that several choices of rric, f^R and fip niay work equally well. In 
fact the cross sections increase by decreasing ^p, fiji or rric, so changes in the three 
variables can be played against each other to obtain practically the same results. We 
present here just one such choice. 

We choose fiR = rrix, fip = 2mT for all sets, and 

rric = 1-185 GeV for MRS Rl, (1) 

rric = 1.31 GeV for MRS R2, (2) 

rric = 1.24 GeV for CTEQ 3M, (3) 

rric = 1-27 GeV for CTEQ 4M. (4) 

We fit rric, fJ'R, and /ij? to the latest available data on charm production |, |^, |10| 
in proton-nucleon and pion-nucleon collisions. We use mainly the data on pN collisions, 
which are more relevant to us, but examine also the nN data to see how well our choice 
of parameters works there. 

The MNR program calculates the total cross section for cc pair production, acc- 
We converted the experimental data on or D~ production a{D^ , D~), or Z)° 



4 



production a{D^ , D'^), or the same cross sections just for xp > {xp is the Feynman 
x), a+{D^, D^) and cr+(D°, into dec values following [|10 . 



The data we used for the 'calibration' of the MNR program are shown in Table |T] 
and Table ^ 0, ^, |T0[ . These tables also present a comparison of experimental data 
on total inclusive D-production cross sections (converted to total cross sections) 
with those calculated with the MNR program. 

For the data of Table |l|, for pN collisions, the conversion is done using 

(Tec = 1.5 X ^ X D-) + a{D\ (5) 

if cross sections are measured for any xp, oi 

cTec- = 1.5 X 2 X i D-) + D^)] , (6) 

if experimental data are given for xp > Q only. The explanation of the factors in Eqs. 
(iD)® is as follows. The | factors convert single D inclusive into DD pair inclusive 
cross sections. The 1.5 factors are required to take into account the production of Ds 



and Ac (which is included in acs) through the ratios [|T0 



"^""'^ -0.2, -#4t-0.3, (7) 



a{D+,DO) ' a{D+,DO) 

(the same relation also for antiparticles). The factor 2 in Eq. (P) converts from xp > 
to all Xp (i.e. it is (Jcc/(^cc{xp > 0) for the pN case). 

In the case of ttN collisions (Table |^) the factor 2 in equation (^ is replaced by 
1.6, which is the value of (Jcc/o'cc{xp > 0) when a pion beam is used. 

Table |l| explains our choice of nic values. The rric values in Eqs.(|l]),(||),(|^) and @ 
reproduce well the central values of the pN charm inclusive total cross sections |^, 
using the program with the four different PDFs. 

In Table ^ we also present a similar analysis for tcN collisions, using only MRS 
Rl for simplicity. In this case slightly higher values of rric fit the vrA^ data |^, |T0| a 
bit better, while rric = 1.185 GeV, the value we take with the MRS Rl PDF, fits the 
pN data [0, H, a bit better. Notice that for the pions we used a different PDF, 
SMR2 the same used in Refs. 0, (obviously not used in our calculations of 



atmospheric fluxes). We present the nN data just for completeness, to show that they 
too are reasonably well fitted with our choice of parameters. These other values of 
in Table 0well reproduce the 7r=^A^ data at 250 GeV [§ and the vr'A^ data at 350 GeV 
(which seem a bit too low with respect to the data at 250 GeV). Even if each value 
of rric reproduces best each total cross section, all three provide reasonable fits to all 
data, as can be seen also in the Figs. 1-3. 

In Figs. 1-3 we present total and differential cross sections calculated with the MNR 
program and compared to the experimental data. As a way of example, we describe 
our fits for MRS Rl only. 

Fig. la shows the fit to pN total cross sections (converted into acc values as described 
above). In addition to the experimental value of Table |l] — which is the fundamental 
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one, since it's the experiment whose differential cross sections we want also to fit — 
we added other experimental points coming from previous experiments (for details see 
[|T0[]). For pN the ttlc = 1.185 GeV is the best choice. 

Fig. lb shows the same for vrA^ collisions. Here, as explained before, values of 
rric = 1.25 GeV or rric = 1.31 GeV are a better choice. Again we added here for 
completeness other experimental points coming from previous experiments . 

Fig. 2ab shows fits to D-inclusive differential cross sections. In this figure the 
theoretically obtained dacc/dxp and dacc/dp^ were converted into D-cross sections, 
with no extra factors. Fig. 2ab presents the data of the E769 collaboration p| for pN 
and ttN at 250 GeV. In these cases the differential cTcc cross sections are converted 
into single inclusive ones (by a factor of 2) and then into cross sections for production 
of and D$ (by a factor of 1.2/1.5, see Eq. (0)) for the E769 data. For 

example, 

|L,fl*fl«.fl«.fl|).Hx2x|;£ (8) 

for Fig. 2a (and similar factors for da/dp\ for Fig. 2b). The fit to the da/dp^ pN 
data in Fig. 2b seems to be a bit too low, but it is not very different from the fit shown 
in Fig. 2 of reference IQ. The predicted da/dp^ are not sensitive to differences in m,c 
that are instead more noticeable in da/dxp- 

Fig. Sab presents the vrA^ data at 350 GeV of the WA92 collaboration |^ in a way 
similar to Fig. 2ab. In these cases the differential cross sections are converted into 
a single inclusive ones (by a factor of 2) and then into cross sections for production of 
D^, D° and D° only (by a factor of 1.0/1.5, see Eq. (0)) for the WA92 data. Similar 
conclusions can be drawn: for pions rric = 1.31 GeV is the best choice in this case. 

We have performed the same analysis with MRS R2, CTEQ 4M and CTEQ 3M, 
even if we do not show here any of the fits. The results for total and differential cross 
sections were similar to those shown for the MRS Rl, the only difference being the 
choice of rric- 

In conclusion, we obtain good fits to all data on charm production with one choice 
of and iTLc for each PDF, without other normalizations. 



2.2 Choice of PDFs 

Consider the collision of a cosmic ray nucleus of energy E per nucleon, with a nucleus 
of the atmosphere in which charm quarks of energy Ef. are produced, which decay into 
leptons of energy Ei (in the lab. frame, namely the atmosphere rest frame). Due to 
the steep decrease with increasing energy of the incoming flux of cosmic rays, only 
the most energetic charm quarks produced count for the final lepton flux, and these c 
quarks come from the interactions of projectile partons carrying a large fraction of the 
incoming nucleon momentum. Thus, the characteristic x of the projectile parton, that 
we call Xi, is large. It is xi ~ O(10~^). We can, then, immediately understand that 
very small parton momentum fractions are needed in our calculation, because typical 
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partonic center of mass energies are close to the cc threshold, 2mc ~ 2 GeV (since 
the differential cross section decreases with increasing s), while the total center of mass 
energy squared is s = 2m nE (with rriN the nucleon mass, m^v — 1 GeV). Calling X2 
the momentum fraction of the target parton (in the nuclei of the atmosphere), then, 
XiX2 = s/s = Aml/{2mNE) ~ GeV/E. Thus, ~ O(GeV/0.1 E), where E is the 
energy per nucleon of the incoming cosmic ray in the lab. frame. The characteristic 
energy E^ of the charm quark and the dominant leptonic energy Ei in the fluxes are 
El c::^ Ec^ O.IE, thus X2 ^ 0(GeV/ Ei), as mentioned above. 

For X > 10"^ {E < 10^ TeV), PDFs are available from global analyses of existing 
data. We use four sets of PDFs. MRS Rl, MRS R2 H and CTEQ 4M incorporate 



most of the latest HERA data and cover the range of parton momentum fractions 
X > 10"^ and momentum transfers > 1.25 - 2.56 GeV^. MRS Rl and MRS R2 
differ only in the value of the strong coupling constant at the Z boson mass: in MRS 
Rl a,(M|) = 0.113, and in MRS R2 a,(M|) = 0.120. The former value is suggested 
by "deep inelastic scattering" experiments, and the latter by LEP measurements. This 
difference leads to different values of the PDF parameters at the reference momentum 
Ql = 1.25 GeV^ where the QCD evolution of the MRS Rl and R2 PDFs is started. 
The CTEQ 4M is the standard choice in the MS scheme in the most recent group of 
PDFs from the CTEQ group (as(M|) =0.116 for CTEQ 4M). We also use an older 
PDF by the CTEQ group, namely the CTEQ 3M |T9|, only for comparisons with [ p!7[| , 
where it is used as the main PDF. 

For X < 10~^ {E > 10^ TeV), we need to extrapolate the available PDFs. For 
X <^ 1, all these PDFs go as 

x/i(x,Q2)^A,a;-"»(«'), (9) 

where i denotes valence quarks u^.d^, sea quarks S, or gluons g . The PDFs we used 
(except the older CTEQ 3M) have \s{.Ql) 7^ ^giQ^)^ contrast to older sets of PDFs 
which assumed an equality. As x decreases the density of gluons grows rapidly. At 
X ~ 0.3 it is comparable to the quark densities but, as x decreases it increasingly 
dominates over the quark densities, which become negligible at x ^ 10^^. 

We need, therefore, to extrapolate the gluon PDFs to x < 10"^. Extrapolations 
based on Regge analysis usually propose X5'(x) ~ x~^ with A ~ 0.08 while evo- 
lution equations used to resum the large logarithms as ln(l/x) mentioned above, such 
as the BFKL (Balitsky, Fadin, Kuraev, Lipatov ||2^) find also xg{x) ~ x~^ but with 



A ~ 0.5 A detailed analysis of the dependence of the neutrino fluxes on the 

low X behavior of the PDFs will be given in another publication [|^. As mentioned 



above, in the present paper our goal is to compare NLO to BORN simulations, for 
which we use a simplified extrapolation at low x of the gluon PDF, which is somewhat 
in between the two extreme theoretical behaviors described above. For MRS R1-R2 
and CTEQ 4M we take a linear extrapolation of lng{x) as a function of Inx, in which 
we took lng{x) = —{^giQ"^) + l)lnx + InA^, where Xg{Q^) was taken as its value at 
X = 10~^, the smallest x for which the PDFs are provided; for the CTEQ 3M we used 
a polynomial approximation which is included in the PDF package. 
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3 Simulation of particle cascades in the atmosphere 



We simulate the charm production process in the atmosphere and the subsequent 
particle cascades, by modifying and combining together two different programs: the 
MNR routines [0 and PYTHIA 6.115 [gg. 



The MNR program was modified to become an event generator for charm produc- 
tion at different heights in the atmosphere and for different energies of the incoming 
primary cosmic rays. 

The charm quarks (and antiquarks) generated by this first stage of the program 
are then fed into a second part which handles quark showering, fragmentation and 
the interactions and decays of the particles down to the final leptons. The cascade 
evolution is therefore followed throughout the atmosphere: the muon and neutrino 
fluxes at sea level are the flnal output of the process. 

In this section we give a brief description of the main parts of the simulation. 
Even if our program is completely different from the one used by TIG, because it is 
constructed around the MNR main routines, nevertheless we keep the same modeling 
of the atmosphere and of the primary cosmic ray flux as in TIG and the same treatment 
of particle interactions and decays in the cascade. 

Our main improvement is the inclusion of a true NLO contribution for charm pro- 
duction (and updated PDFs), so we keep all other assumptions of the TIG model in 
order to make our results comparable to those of TIG. We study the effect of modifying 



some of their other assumptions elsewhere 1 24 



3.1 The model for the atmosphere 

We assume a simple isothermal model for the atmosphere. Its density at vertical height 
h is 

p{h) = ^e-^'''\ (10) 
ho 

where the scale height = 6.4 km and the column density Xq = 1300 g/cm^ at 
h = are chosen as in TIG, to flt the actual density in the range 3 km < h < 40 km, 
important for cosmic ray interactions. Along the vertical direction, the amount of 
atmosphere traversed by a particle, the depth X, is related to the height h simply by 

X= / p(h')dh' = Xoe-''/^°. (11) 
Jh 

The atmospheric composition at the important heights is approximately constant: 
78.4% nitrogen, 21.1% oxygen and 0.5% argon with average atomic number (A) = 
14.5. 
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3.2 The primary cosmic ray flux 



Following TIG W^, we neglect the detailed cosmic ray composition and consider all 



primaries to be nucleons with energy spectrum 
nucleons 



(E,0) 



cm^ s sr GeV /A 



f 1.7(E/GeV)-2-7 for ^ < 5 10^ GeV , . 
\ 174(E/GeV)-3 for ^ > 5 10^ GeV ^ ' 



The primary flux is attenuated as it penetrates into the atmosphere by collisions against 
the air nuclei. An approximate expression for the intensity of the primary flux at a 
depth X is (see again) 

ME,X) = e-''/''-^^^ ME,0) . (13) 

The nuclear attenuation length Atv, defined as 

has a mild energy dependence through Znn and \n- Here Z^^ is the spectrum- 
weighted moment for nucleon regeneration in nucleon-nucleon collisions, for which we 



use the values in Fig. 4 of Ref . |TT| . And A at is the interaction thickness 



XN{E,h) = - ^1^) 

22A(^NA{E)nA{h) 

where nA{h) is the number density of air nuclei of atomic weight A at height h and 
a]siA{E) is the total inelastic cross section for collisions of a nucleon N with a nucleus 

This cross section scales essentially as A^/^, since for the large nucleon-nucleon 
cross sections we deal with, the projectiles do not penetrate the nucleus. So we set 

{E) = A^/^aNN{E). For aNN{E) we use the fit to the available data in Ref. [^^ 



Using our height independent atmospheric composition, we simplify Eq. ([T5|) as follows, 

A^(E, h) = ^ = 2.44 — ^ . (16) 

Here ( ) denotes average and u is the atomic mass unit, that we write as 

u = 1660.54 mb g/cm^. (17) 

We therefore find that in our approximations XnIE) is independent of height. 

*We recall that the elastic cross section contributes negligibly to the primary flux attenuation 
because the average elastic energy loss is very small, less than 1 GeV at the high energies we consider. 
This can be seen using the differential elastic cross section daei/dQ^ — {daei/dQ^)Q2^Qexp{—bQ^) 
with b — [7.9 + 0.9 lnp;afc]GeV^^, with piab in GeV Here Q is the momentum transfer of the 

colliding proton of incoming momentum piat and mass M . The mean energy loss is the mean value of 
Q^/2M (here M is the target proton mass) namely {l/2Mb) = 67MeV/(l + 0.1 ln(p,ab/GeV)). This 
is 46 MeV at £' = lOOGeV, and smaller at higher energies. 
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3.3 Charm production with MNR routines 

As we remarked before, the modified MNR routines are the first stage of our simulation. 
For a given energy E of a primary incoming proton in the lab system, i.e. in the 
atmosphere reference frame, we generate a collision with a nuclear target at rest in 
the atmosphere, activating the MNR routines (primary event, pN collision, with = 
{p + n)/2) . 

These routines generate total and differential cross sections through a VEGAS 
integration, which creates a large number of 'subevents', each one with a particular 
weight, which in the original MNR program are summed together to calculate the final 
cross sections. 

It is easy to modify the program so that each of these subevents (together with 
its weight) can represent the production of a charm c (or of a cc pair, or cc gluon, 
etc.) with given kinematics in any particular reference frame of interest. The original 
MNR routines can calculate single differential cross sections, in which the kinematics 
of only one final c quark is available, and double differential cross sections, in which the 
full kinematics of the cc pair (plus an additional parton in NLO processes) becomes 
available, for each subprocess. We have used both these possibilities. We will refer 
to them as 'single' and 'double' modes. The 'single' is the mode we use to obtain all 
our results. We use the 'double' mode only to compare the results of the independent 
fragmentation model used in the evolution of cascades in the 'single' mode, with the 
more reliable string fragmentation model, which can only be used in the 'double' mode, 
as we explain below. 



The MNR program |T5|, |T6[ contains all BORN and NLO processes. In the 'single' 
mode we can generate the following processes, with only the kinematics of the c quark 
available, 
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cX; qq cX (BORN) gg cX; qq cX; qg cX (NLO) 



where q represents any light quark or antiquark. In the 'double' mode we have the 
following processes 

gg cc; qq — » cc (BORN) gg — > ccg; qq ccg; qg ccq (NLO) (19) 

for which the kinematics of all the outgoing partons is fully determined for each 
'subevent'. 

All the kinematical variables of the partons in the final state constitute the input 
for the next stage of the program, described in the next subsection. 

An important characteristic of the first stage is that, besides rric, fiR, and fip, we 
can select any desired PDF to be used with the charm production routines. We have 
updated the set of PDFs in the original MNR program. 

According to the discussion of Sect. |, we use the MRS Rl, MRS R2, CTEQ 3M 
and CTEQ 4M parton distribution functions, together with the values of mc, fiR, and 
I^F in Eqs. (^). 
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As a concrete example of the integrals performed in our program, here we write the 
differential flux 0^ of muons (namely of /i^ + with energy {fi stands here for 
li'^ or in the 'single' mode (0^ has units cm~^ s~^ sr^^ GeV 



/•OO /"OO 

,{E^) = dE dhME,X{h))J2nAih) 

•lEu Jo . 



X 



dEr. 



+ 



da{pA cY; E, E^ 
dE, 



J MNR I 



dNf,{c /i; Ec, E^, h) 
dE„, 



PYTHIA 

(20) 



Here nA{h) is the number density of nuclei of atomic number A in the atmosphere, E 
is the energy of the primary cosmic ray proton, E, the energy of the charm produced 
in the collision pA cY {Y here stands for anything else). Using the relation 
da{pA cY)/dEc = A da{pN cY)/dEc, the sum over A becomes J^a'i^a^A = 
p{h)/u. Using dX = —p{h)dh, Eq. (p!3|), and normalizing to one the distribution in 
depth X, (pa becomes 



,{E,) = dE dX ME,X=0) ——— 



f{h)AN{E) 



u 



(21) 



where, from Eqs.(|T|) and (|T6D, Ajy/u = 2.44[oriVAr(l - Znn)Y^ and 



f{h) 



dEr 



da{pN cY; E, E, 
dEr 



J MNR. 



dN^{c fi]Ec,Ef„h) 
dE„ 



(22) 



J PYTHIA 



Here the factor of 2 accounts for the muons produced by c (only c quarks are used in the 
program for simplicity); the pN inclusive charm production cross section is computed 
with the MNR program (here are the integrations over the PDFs and partonic cross 
sections) and the last square bracket is the number of muons of energy E^ which reach 
sea level, produced in the cascades simulated by PYTHIA. Each cascade is initiated by 
a c quark (in the 'single' case) of energy E, and momentum k (provided by the MNR 
routines) at a height h chosen through a random number R homogeneously distributed 
between and 1, which gives the value of the X probability distribution in Eq. (|2T|), 
namely R = e'^^^^^^K 



3.4 Cascade evolution with PYTHIA routines 

The parton c (or partons in the 'double' case) generated by the first stage, namely 
by the MNR routines, are entered in the event list of PYTHIA and they become the 
starting point of the cascade generation. 

PYTHIA first fragments the c quark (in the 'single' mode, or all the partons in the 
'double' mode) into hadrons, after showering, which can be optionally shut off. The 
charm quarks hadronize into D^, D^, Df and Ac. We used here the Peterson 
fragmentation function option. For each hadron produced, a simple routine added to 
PYTHIA decides if the hadron interacts in the atmosphere (loosing some energy) or 
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decays. This is the same approach as in TIG. PYTHIA follows in this way the cascade 
in the atmosphere and populates the histograms of muons and neutrinos as a function 
of their different energies. We mention here a few important technical details. The 
'single' and 'double' modes described before use different fragmentation models. In the 
'single' mode only one c quark is available and is entered at the beginning of the event 
list (with its energy and momentum in the partonic CM reference frame). In this case 



PYTHIA uses the 'independent fragmentation' model (see p5| for details). We only 
include c quarks and at the end multiply the result by a factor of two to account for 
initial c quarks. 

In the 'double' mode, instead, which we only use at the LO, we start with two (cc) 
partons in the event list. In this case we opt to use the 'string fragmentation' model 
(Lund model, |^). This model generally gives better results than the independent 
fragmentation, in which energy and momentum conservation have to be imposed a 
posteriori and whose results depend on the reference frame used, which empirically is 
chosen to be the partonic CM frame. To impose energy and momentum conservation 
in the independent fragmentation, we used the option (MSTJ(3)=1, see again [^) in 
which particles share momentum imbalance compensation according to their energy 
(roughly equivalent to boosting events to CM frame) but we have convinced ourselves 
that the results do not depend much on the way of imposing energy and/or momentum 
conservation, because trial runs with different options have given similar results for the 
fluxes. 

Even if independent fragmentation is in general less desirable than string fragmen- 
tation, we use the 'single' mode as our main choice. The main reason to use the 'single' 
mode is that the simulations run in acceptably short times (a few days) on the SUN 
computers we use, while giving results practically identical to the 'double' mode in 
the comparisons we have made (see Fig. 6c). The simulation of the cascades in the 
'double' mode takes between five and ten times longer. We tested the goodness of 
the independent fragmentation by comparing the outcome of fluxes computed at the 
Born level, in which the charm fluxes at production are identical (we put one c in the 
atmosphere and multiply the outcome by two to account for the c in one case, and we 
put cc in the atmosphere, instead, in the second case) and the sole difference in both 
modes is due to the different fragmentation models used. The results were extremely 
close (at Born level the difference is less than 5%, at energies above lO^GeV), as can 
be seen in Fig. 6c. 

Apart from the mentioned differences between the 'single' and 'double' modes, the 
simulations then proceed basically in the same way in both modes. For each of the 
'subevents', i.e. for each set of initial parton(s) put in the event list, a certain height 
in the atmosphere is randomly chosen as explained above, this being the position at 
which the partons are generated from the initial proton-nucleon collision. This random 



height h is generated in a way similar to TIG (see Ref. ^^), but different, because we 
include a correction for nucleon regeneration in nucleon-nucleon collisions by using A^r, 
the nuclear attenuation length, in Eq. (^) instead of Aat , the interaction thickness 
(see Eqs. (p!4D , ([l5|) and (|l^)).The only difference compared to TIG (see Eq. (15) in 
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the last paper of Ref. |TT|) is the inclusion of the (1 — Z^n) correction term. This 
was done because we could not include regenerated protons directly in our simulation 
of the cascades, since events and subevents are now created by the MNR routines and 
not by PYTHIA, as it was in TIG. 

When parton showering is included at the beginning of the cascade simulation 
performed by PYTHIA, some double counting is present. The double counting appears 
when a LO diagram, for example gg — > cc, with a subsequent splitting contained in 
PYTHIA, for example c gc is summed to NLO diagram, gg gcc with the same 
topology, as if both diagram were independent, when actually the NLO contains the 
first contribution when the intermediate c quark on mass shell. We have not tried to 
correct this double counting but have instead confronted the results obtained including 
showering (our standard option) with those excluding showering (in which case there 
is no double counting) and found very similar leptonic fluxes (see Fig. 6b). 

The particles generated after the initial hadronization are then followed throughout 
the atmosphere and PYTHIA evolves the cascade with the same treatment of inter- 
actions and decays proposed by TIG. The final number of muons and neutrinos at 
sea level is therefore calculated considering all the 'subevents', each with its respective 
weight Wi from the MNR program, which produce the final particles through all the 
possible decay channels of charmed particles decaying into prompt leptons. Since only 
the decay modes of charmed hadrons going into /i or z/^ or i^e are left open in the sim- 
ulation, and there are essentially just 2 modes for each charmed particle (for example: 
Ve + anything , with branching ratio = 0.172; — > /i"*" i/^ + anything, 
with branching ratio = 0.172; all other channels closed), the branching ratios for each 
of these modes is fictitiously taken by PYTHIA to be 1/2 and need to be normalized by 
multiplying by the actual branching ratio (0.172 for the example above) and dividing 
by 1/2. Besides, since not all events are accepted by PYTHIA to generate a complete 
cascade, the result is normalized by dividing by the sum of all the weights of accepted 
events and multiplying it by the total c inclusive cross section. 

3.5 Summary 

To summarize, our computation of the final fluxes is organized as follows. 

• An external loop over the primary energy E generates an integration over E in 
the range 10^ - lO^GeV". 

• For each primary energy E, the MNR routines generate 'subevents' with weight 
Wj, for all the LO and NLO processes. 

• Each subevent is assigned a random height (so that implicitly an integration over 
h is performed) and all this is passed to PYTHIA as a definite set of parton(s) to be 
put at the beginning of the event list. 

• For each of these 'subevents', PYTHIA treats showering (in our standard option), 
hadronization and evolution of the cascade in the atmosphere, and generates the final 
leptons. 
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• For each decay channel of interest, the produced leptons are weighted with Wi 
and then summed into the final fiuxes. 

4 Neutrino and muon fluxes 

Figs. 4-6 show the results of our simulations. Fig. 4 shows the total inclusive charm- 
anticharm production cross sections Ucc, and the K factor for c production, namely the 
ratio between the NLO and Born cross sections, Kc = c^'^'^/cr^"^"', for the four PDFs 
we consider and for TIG. Fig. 5 shows our main results obtained with our default choice 
of options: a 'single' mode calculation including the contributions from all processes 
in Eq. ([TsD and with parton showering included in the cascade simulation performed 
by PYTHIA). Finally Fig. 6 shows the relative importance of the processes included 
in the fiuxes and a comparison of the 'single' and 'double' modes and of the 'on' and 
'off' showering options. 

In Fig. 4a, the total inclusive charm-anticharm production cross sections cTcc are 
plotted over the energy range needed by our program, E < 10^^ GeV, for our four 
different PDFs. They were calculated using the MNR program, with the 'calibration' 
described in Sect. ^, up to the NLO contribution. For comparison, we also show the 
cross section used by TIG and the Born (LO) contribution for one of the PDFs, MRS 
Rl. We see in the figure that all our cross sections agree at low energies, as expected 
due to our 'calibration' at 250 GeV, and are very similar for energies up to 10^ — 10'' 
GeV. At higher energies they diverge, differing by at most 50% at the highest energy 
we use, 10^^ GeV. In fact, at energies beyond 10'' GeV, the CTEQ 3M cross section 
becomes progressively larger than the CTEQ 4M and MRS R2 cross sections, which 
are very close to each other. The MRS Rl becomes on the contrary progressively lower 
than the other three. 

We see in Fig. 4a that for energies above 10^ GeV our cross sections are considerably 
higher than the one used by TIG. This difference can be traced in part to the use by 
TIG of an option of PYTHIA by which the gluon PDF is extrapolated to x < 10~^ 
with A = 0.08, while all the PDFs we use have a higher value of A ~ 0.2 — 0.3. And 
in part to TIG scaling the LO cross sections obtained with PYTHIA by a constant 
K factor of 2, while at large energies the K factor is actually larger than 2 by about 
10-15% (see Fig. 4b). 

In Fig. 4b we explicitly show the K factor for c production, namely the ratio between 
the NLO and Born cross sections, = for our PDFs and for TIG. All 

the Kc values are around the usually cited value of 2 for most of the intermediate 
energies, but are larger at the lowest energies and also at the highest energies (except 
for CTEQ 3M), and they all are within about 15% of each other. 

Fig. 5 contains three sets of figures, one for each lepton: /i, z/^ and Ug. The left figure 
of each set shows the £'^-weighted vertical prompt fiuxes, for all our PDFs up to NLO 
(labelled 'NLO') and, as an example, the LO (labelled 'BORN') for MRS Rl, together 
with the total fiuxes up to NLO of TIG, both from prompt and conventional sources 
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(dotted lines). The right part of each set shows the corresponding Ki value (where 
/ = /i, v^, Ve)i i-e. the ratio of the total NLO flux to the Born flux of the figure on the 
left. The figures show that our fluxes are higher than those of TIG for E > 10^ GeV. 
Leaving apart differences in the two simulations that cannot be easily quantified, this 
discrepancy can largely be explained by the different cross sections used by TIG and 
us: the TIG cross section is lower than ours for E > 10'^ GeV. Using a value of A 
similar to TIG (A ~ 0) at small x, we obtain fluxes similar to those of TIG at energies 
above 10^ GeV |2|. 

In particular, our fluxes are all larger than TIG by factors of 3 to 10 at the highest 



energies, what puts our fluxes in the bulk-part of previous estimates (see Refs. [|12 



T3| , |T4|, ^). There is an evident dependence of the fluxes on the choice of PDF. It is 
remarkable that MRS R2 and CTEQ 4M give very similar results. Those of the MRS 
Rl become lower and those of the older CTEQ 3M PDF become higher as the energy 
increases (both differing by about 30-50% at the highest energies with respect to the 
MRS R2-CTEQ 4M fluxes). This is due to the intrinsic differences of the PDF packages 
used and the consequent different extrapolated values of A at small x or high energies. 
The CTEQ 3M fluxes were included to compare our results with those of Ref. 



1^. We find our CTEQ 3M results to be close to those of Ref. |T^, in spite of 



the very different approaches used in the two calculations. Our fluxes lie between 
the two curves for CTEQ 3M shown in Fig. 8 of Ref. Jl^, corresponding to different 



choices of renormalization and factorization scales. Our fluxes are lower (by 30-40% 
at lO^GeV), than the main CTEQ 3M choice of Ref. fl^ (solid line of their Fig.8), 
which is calculated using values of fiR, fip and rric similar to ours. Our cross section for 
charm production, for the CTEQ 3M case, is essentially equal to the one used in Ref. 
[0 (shown in their Fig. 2), so the discrepancies in the final fluxes are to be explained 
in terms of the differences in the cascade treatment. It is very difficult to trace the 
reasons for these differences. 

We also see in the figures that, for each PDF, the fluxes for the different leptons are 
very similar: those for neutrinos and Ue are essentially the same, those for muons are 
only slightly lower (about 10% less at the energies of interest). Also the Ki's don't differ 
much for the three leptons, apart from some unphysical fluctuations especially evident 
at the highest energies. Even if they differ the PDF, they all show a similar energy 
dependence, namely they increase at low energies and sometimes at high energies also. 
This behavior is also similar to that of the Kc factors in Fig. 4b, but with a weaker 
overall energy dependence, as expected, since the leptons of a given energy result from 
c quarks with a range of higher energies. 

The Ki factors are all within the range 2.1 — 2.5: they are approximately 2.2 for 
MRS Rl, 2.4 for MRS R2 and CTEQ 4M, and 2.3 for CTEQ 3M. Thus, our analysis 
shows that evaluating the lepton fluxes only at the Born level, and multiplying them 
by an overall Ki factor of about 2.2 — 2.4 (i.e. 10 to 20% larger than the value of 2 
used by TIG|^, can be good enough to evaluate the NLO fluxes within about 10%. 
Thus we find the approach used by TIG, who multiplied the LO fluxes obtained with 



t We note that in the original TIG model there is no distinction between Kc and Ki factors since 
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PYTHIA by two, essentially correct, except for their relatively low K factor and the 
discrepancies existing even at Born level between our fluxes and those of TIG. In fact, 
as we mentioned previously, the differences between our final results and those of TIG 
depend mostly on the different total inclusive c cross sections, which can be traced to 
the extrapolation of the gluon PDF at small x rather than to the K factor. Possible 
causes of the different results due to the intrinsic differences of the computer simulations 
cannot be easily quantified. 

In Fig. 6 we address three issues. First, we show that the fluxes can be obtained 
within about 30% with just the gluon-gluon process. This would speed up the simu- 
lations and, when using the MNR program, would give (contrary to intuition) higher 
fluxes than those actually derived from all processes. Secondly, we show that the fluxes 
obtained including or excluding showering in the simulation made by PYTHIA (we in- 
cluded showering in our standard options) do not differ significantly. The third issue 
we deal with is the difference between the 'single' and 'double' modes described in 
Sect. 3. We show that at LO the results from a 'double' mode calculation coincide 
with those of the much shorter 'single' mode, that we use in all our calculations. Let 
us deal with these three issues in turn. 

In Fig. 6a we show, for a given PDF, the MRS Rl, the relative importance of 
the different processes contributing to the final fluxes. The solid line is the total flux 
obtained as the sum of all the processes of Eq. ([T8|) and the dotted line shows the result 
of only gluon-gluon fusion {gg)., the sum of Born [gg) and pure NLO (excluding Born) 
gg processes. Also shown are the separate contributions only at the Born and at the 
NLO (excluding LO) of both gg and quark- ant iquark (gg) fusion, what clearly shows 
that gg dominates. This is to be expected because the gluon PDF is either much larger 
than (for x < 0.1) or comparable to (for x ~ (9(0.1)) the quark PDFs. The figure 
plots the absolute value of the quark-gluon {qg) terms because, for the values of the 
factorization scale that we employ in our calculations, these terms are negative. This 
is due to the way the original MNR calculation is subdivided into processes. In fact, 
in the MNR program, a part of the quark-gluon contribution to the cross sections is 
already contained in other processes, and must be subtracted in the processes labelled 
as qg. The amount subtracted depends on the factorization scale and may drive 
the qg contribution negative. Roughly speaking, if /x^? is small the qg term is positive, 
otherwise (as in our case) the term is negative. The absolute value of the qg term is 
in between the gg and the gg terms, what makes negative the sum of all the processes 
different from gg. Thus, gluon-gluon processes alone give a result slightly larger than 
the total, by about 30%. 

In Fig. 6b we check the effect of shutting off the showering option available in 
PYTHIA. We study only one specific case, the MRS Rl. The overall effect is minimal: 
the exclusion of showering slightly increases the energy of the parent charmed hadrons 
and therefore causes the final fluxes of lepton daughters to move towards higher en- 
ergies; the effect is barely noticeable and just slightly more important for the Born 

only the Born level is considered. Their K — 2 factor is just a multiplicative constant which can be 
considered either a. oi a, Ki. 
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fluxes (the overall difference is about 5%). When showering is included some double 
counting occurs, whose effect must be smaller than the difference between the results 
with showering on and off (since in this case no double counting occurs). 

Finally in Fig. 6c we confront the 'single' and 'double' modes of the program, for 
just one PDF, MRS Rl, at Born level. At this level, the calculation of the charm 
flux at production is identical (we obtain the fluxes from c and multiply by two at 
the end to account for the c in one case, and we obtain the fluxes directly from cc 
in the other). So, what is actually compared in the two modes at the Born level is 
the fragmentation model: independent fragmentation in the 'single' mode and string 
(Lund) fragmentation in the 'double' mode. The results from both modes at the Born 
level are almost identical: as already remarked the difference is less than 5% for energies 
above 10^ GeV. 

5 Conclusions 

We have used the actual next-to-leading order perturbative QCD calculations of charm 
production cross sections, together with a full simulation of the atmospheric cascades, 
to obtain the vertical prompt fluxes of neutrinos and muons. 

Our treatment is similar to the one used by TIG, except for the very important 
difference of including the true NLO contribution, while TIG used the LO charm pro- 
duction cross section multiplied by a constant K factor of 2 to bring it in line with the 
next-to-leading order values. The main goal of this paper is to examine the validity of 
TlG's procedure by computing the ratio of the fluxes obtained with the NLO charm 
production cross section versus those obtained with the LO cross section. 

These ratios, the Ki factors are between 2.1 and 2.5 for the different gluon PDFs 
in the energy range from 10^ to 10^ GeV (see Fig. 5). Consequently, our analysis 
shows that evaluating the lepton fluxes only at the Born level, and multiplying them 
by an overall factor of about 2.2 — 2.4, slightly dependent on the PDF, can be good 
enough to evaluate the NLO fluxes within about 10%. Therefore, we find the approach 
used by TIG (i.e. multiplying the LO fluxes by two) essentially correct, except for their 
relatively low K factor. We find different lepton fluxes than TIG, but this is mostly due 
to the discrepancies, even at Born level, between our charm production cross sections 
and TIG'S. 

In fact, the prompt neutrino and muon fluxes found by TIG were lower than the 
lowest previous prediction. We find here instead fluxes in the bulk part of those previous 
predictions. This difference can be traced largely to the use by TIG of an option of 
PYTHIA by which the gluon PDF is extrapolated for x < 10~^ with A = 0.08, while 
all the PDFs in this paper have a higher value of A ~ 0.2 — 0.3. Using a value of A 
similar to TIG (A ^ 0) we obtain fluxes similar to those of TIG, at energies above 10^ 
GeV m- 
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FIGURE CAPTIONS 



Fig. 1 Comparison of experimental data for Ucg with MNR predictions for different rric 
values: (a) in pN collisions (flOl? Table 0), (b) in nN collisions {^^^, Table |^) 



(PDF: MRS Rl). 

Fig. 2 Comparison of differential cross sections for {D~^, D~, D", D^, Dg and Dg) pro- 
duction, calculated using MNR at different rric values, with E769 data for pN and 
vrA^ [|]: (a) da/dxp, (b) da/dp'^ {xp > 0) (PDF: MRS Rl). 

Fig. 3 Comparison of differential cross sections for {D^ , D~ , D'^ , D^) production, cal- 
culated using MNR at different rric values, with WA92 data for vrA^ [^: (a) 
da/dxp, (b) da/dp^ {xp > 0) (PDF: MRS Rl). 

Fig. 4 (a) Total cross sections for charm production dec up to NLO, for different PDFs, 
compared to the one used in the TIG model |]ll| (for MRS Rl we also show the 
Born cross section), (b) Related Kc factors. 

Fig. 5 i?^-weighted vertical prompt fluxes, for different PDFs, at NLO (for MRS Rl 
we also show the Born flux), for the three types of leptons considered, compared 
to the TIG |Tl| conventional and prompt fluxes (left figures) and the related Ki 



factors for each case (right figures). 

Fig. 6 (a) Contributions of the different Born and NLO processes to the total E^- 
weighted vertical prompt fluxes, (b) Comparison of the fluxes with or without 
the showering option, at Born and NLO. (c) Comparison of the fluxes calculated 
in the 'single' or 'double' mode, at Born only (PDF: MRS Rl). 
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1 Q c _|_ 9 9 


iO.041: 

mc = 1.185 GeV 


iVirto SXL 


pN 

E769 


250 


(T+(D+,L)-) = 3.3 ±0.4 ±0.3 
CT+(L>°,D°) = 5.7 ± 1.3 ±0.5 


13.5 ±2.2 


13.43 
= 1.31 GeV 


MRS R2 


pN 

E769 


250 


(T+(D+,D^) = 3.3 ±0.4 ±0.3 
a+{D°,D^) = 5.7 ± 1.3 ±0.5 


13.5 ±2.2 


13.59 
= 1.27 GeV 


CTEQ4M 


E769 


250 


a+{D+,D-) = 3.3 ±0.4 ±0.3 
a+{D^,D^) = 5.7 ± 1.3 ±0.5 


13.5 ±2.2 


13.45 
rric = 1.24 GeV 


CTEQ3M 



Table 1: Data on total cross sections for charm production for pN collisions, from E769 
experiment, have been converted to cc cross sections and compared to the predictions 
of the MNR program running at slightly different values of the charm mass tcIc, using 
different PDFs. 





Beam 
Energy 
(GeV) 


cr+(a;F > 0) 


acc-(EXR) 


acc-(MNR) 
(/ife) 

mc = 1.185GeV 


^cc-(MNR) 

(/ifo) 
= 1.250GeV 


TT-N 
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D+,D- : 1.7 ±0.3 ±0.1 
: 6.4 ±0.9 ±0.3 


9.7 ± 1.2 


14.08 


10.64 


TT-N 

E769 
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D+,D- : 3.6 ±0.2 ±0.2 
: 8.2 ±0.7 ±0.5 


14.2 ± 1.1 


16.54 


12.56 


1T+N 

E769 


250 


D+,D- : 2.6 ±0.3 ±0.2 
D°,D^ : 5.7 ±0.8 ±0.4 


10.0 ± 1.2 


16.54 


12.56 


E769 


250 


D+,D~ : 3.2 ±0.2 ±0.2 
: 7.2 ±0.5 ±0.4 


12.5 ±0.8 


16.54 


12.56 


TT-N 

WA92 


350 


D+,D- : 3.28 ±0.08 ±0.29 
: 7.78 ±0.14 ±0.52 


13.3 ±0.7 


22.22 


17.06 
(13.5 for rric 
=1.31GeV) 



Table 2: Data on total cross sections for charm production for irN collisions, from E769 
and WA92 experiments, have been converted to cc cross sections and compared to the 
predictions of the MNR program running at slightly different values of the charm mass 
mc, using MRS Rl. 
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